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time version of the relation is implemented in Monte Carlo simulations and then used 
to study the aging regime of the ferromagnetic two-dimensional Ising-Glauber model 
quenched from the paramagnetic phase to the ferromagnetic one. Our approach has the 
originality to give direct access to the response function and the fluctuation-dissipation 
ratio. 
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1. Introduction 

The knowledge about out-of-equilibrium processes is far from being as advanced as 
for systems at thermodynamical equilibrium. In particular, the fluctuation-dissipation 
theorem (FDT) which holds at equilibrium is known to be violated out-of-equilibrium. 
This theorem states that at equilibrium the response function R eq (t — s) at time t to an 
infinitesimal field applied to the system at time s < t is related to the time-derivative 
of the two-time autocorrelation function C eq (t — s): 

R*(t-s) = /3£.Cr*(t-s). (1) 

OS 

In the Ising case, the response function reads Rji(t,s) = and the correlation 

function Cji(t,s) = (<x,(t)<7j(s)). Based on a mean-field study of spin-glasses, 
Cugliandolo et al. I\ have conjectured that for asymptotically large times the FDT 
can be generalised by adding a multiplicative factor X(t, s) which moreover depends on 
time only through the correlation function: 

R(t,s) t ~^pX(C(t,s))-?- s C(t,s). (2) 

The quantity /3 eS (t, s) = (3X(C(t, s)) is interpreted as an effective inverse temperature. 
Exact results have been obtained for the ferromagnetic Ising chain 00] that confirm this 
conjecture. Unfortunately, the response function is rarely so easily accessible for more 
complex systems. Both numerically and experimentally, only the integrated response 
function is usually measured by applying a finite magnetic field to the system during 
a finite time. In the so-called TRM scheme, the magnetic field is applied between the 
times and s and the magnetisation is measured at time t. Assuming the equation (J2J) 
valid for any times t and s, one can relate the integrated response function to the 
fluctuation-dissipation (FD) ratio: 

rs rC(t,s) 

x (t,s)= R(t,u)du~p X(C)dC. (3) 

Jo Jc(t,o) 

The FD ratio X(t, s) can thus be obtained as the slope of the integrated response 

function x(t, s) when plotted versus the correlation function C(t,s). This method has 

been applied to the numerical study of many systems: 2d and 3d-Ising ferromagnets [3], 

3d Edwards- Anderson model 0131, 3d and 4d-Gaussian Ising spin-glasses jHj, 2d Ising 

ferromagnet with dipolar interactions [Zj, Heisenberg anti-ferromagnet on the Kagome 

lattice [H] . . . The conjecture (J2J) has also recently been checked experimentally for 

a spin-glass 9 . More details may be found in the reviews ^01 ITTj . However, the 

integrated response function depends linearly on the FD ratio only if the conjecture (0) 

holds, which has not been demonstrated for any of the previously cited systems. We 

will see in the case of the homogeneous Ising model that this approach may lead to 

misinterpretations and erroneous values of X(t, s). The generalisation of the equilibrium 

FDT has recently become an increasingly popular issue. Let us mention two of them: an 

approximate generalisation of the FDT to metastable systems [TJ] (limited to dynamics 

having a transition rate W with only one negative eigenvalue) that has been successfully 
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compared to numerical data for the 2D-Ising model and a generalisation of the FDT for 
trap models [T5| . 

In the present work, we study the dynamics of an Ising-Glauber-like model. In 
the section 1, we describe the model and its dynamics which are studied analytically 
in the section 2. The response function to an infinitesimal magnetic field is exactly 
calculated far-from-equilibrium. It turns out that the response function is no more 
related to a time-derivative of the spin-spin correlation function but to time-derivatives 
of a more complicated correlation function. The equilibrium limit is shown to have the 
usual form. In the section 3, a discrete-time version of this expression is implemented 
in Monte Carlo simulations. Our approach presents several advantages: (i) we can 
compute directly the response function and not only the integrated response function, 
(ii) we obtain the response function to an infinitesimal magnetic field so that we avoid 
non-linear effects due to the use of a finite magnetic field, (iii) the FD ratio can be 
computed without resorting to Cugliandolo conjecture (J2J) and (iv) we can calculate 
the response function R(t, s) and the FD ratio X(t, s) for any time t and s < t during 
one single Monte Carlo simulation. We performed Monte Carlo simulations of the two- 
dimensional homogeneous Ising model quenched at and below the critical temperature 
T c . In both cases, the expected scaling behaviour of the response function in the aging 
regime is well reproduced by the numerical data. The value of the exponent a, still 
controversial, is estimated and the FD ratio is computed. Our estimate of Xoo at T c 
turns out to be compatible with previous work and the scaling behaviour of X(t, s) 
below T c is well reproduced. In both cases, the FD ratio depend on time not only 
through the correlation function. 



2. Our Ising-Glauber-like model 

2.1. Useful relations on Markov processes 

We consider a classical Ising model whose degrees of freedom are N scalar variables 
(jj = ±1 located at the nodes of a <i-dimensional lattice. Let us denote by p({a},t) 
the probability to observe the system in the state {a} at time t. We first define a 
discrete-time Markov chain by the master equation 

p({a}, t + At) = (1 - At)p({a}, t) + At £ W({a'} - {a}, t)p({a'}, t) (4) 

{a'} 

where ^({cr} — > {o~'},t) is the transition rate per unit time from the state {cr} to 
the state {a'} at time t. The condition J2{a'} W({o~} — > {(/},£) = 1 ensures the 
normation of the probability p({a},t) at any time t. The system is not forced to 
make a transition at each time step, i.e. the transition rate may have non-zero diagonal 
elements M^({cr} — > {a},t). In the continuous-time limit At — > 0, the master equation 
(HJ) goes to 

(l + |) KM, t) = E W({a'} {a}, t)p({a'}, t). (5) 
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It is easily shown that the conditional probability, p({a}, t\{a'}, s) with s < t, defined 
by the Bayes relation 

p(W},t) = Y.P(W}'t\W'},s)p({a'},s) (6) 
{*>} 

satisfies the same master equation (J1J too: 

p({a}, t + At\{a'}, s) = (1 - At)p({a}, t\{a'}, s) 

+ At £ ^({a"} - {a},t)p({a"},t\{a'},s) (7) 

or in the continuous-time limit 

(l + |) p({a},t\{a>},s) = W({a"} - {a},t)p({a"},t\{a'}, s). (8) 

Moreover, one can work out a master equation for the time s. It reads 
p{{a},t\{a'},s + At) = (1 + At)p({a},t\{a'},s) 

- At £ w({</} - K},sMM,t|K}, s ) (9) 

and in the continuous-time limit 

'l " p({a},t\{a'},s) =J2W({a f } -> {</'}, s)p({cr}, t|{</'}> s).(10) 



This last equation might be obtained simply by using for example the identity 

£p(M,t) = o. 

When the transition rates do not depend on time, the conditional probability 
p({a}, t\{a'}, s) is a function of t — s only This can be shown easily by introducing the 
matrix notation p({a}, t\{cr'}, s) = ({c}\ p(t,s) \{o~'}). The master equation (JHJ) reads 
then: 

lp(t,s) = (W-l)p(t,s) (11) 
where ({cr}\ W {{&'}) = W({a} — > {er'}). This equation admits the formal solution 

p({a} 7 t\{a'},s)= J2 (W}\e J >~ i)dt ' p({a"},s\{a'},s) 
W"} 

= <M| e^- 1 ^') |{(/}> (12) 

where the initial condition p({a"}, s\{a'}, s) = o~{a"},{a'} has been used. This dependence 
only on t — s, even far-from-equilibrium, will be used latter in the calculation of the 
response function. 

2.2. The model and its dynamics 

The Ising model is defined by its equilibrium probability distribution p eq ({a}) which 
reads with general exchange couplings: 

fWM) = le-WM) = I/£m<^<™. (13) 
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where ferromagnetic couplings correspond to Jki > 0. The condition of stationarity 
({°"}) = leads according to the master equation to a constrain on the 



dt 

transition rates: 

E [Po q (W})W({a'} -> {*},*) - p cq ({a})W({a} -> {</},*)] = 0. (14) 

The equation (JT4*j) is satisfied when the detailed balance holds: 

P^W})W({a'} - W,t) = MMMM -> {</},*). (15) 
This last unnecessary but sufficient condition is fulfilled by the heat-bath single-spin flip 
dynamics defined by the following transition rates: 

1 N 

W({a} - {a'}, t) = -Yl W k ({a} - {a'}) (16) 

k=l 

where the transition rate for a single spin-flip is 
W k ({a} -> {a'}) ^ 



In this last expression, only the single-spin flip — > er^ is allowed. The product of 
Kronecker deltas ensures that all other spins are not modified during the transition. 
After the transition, the spin o~k takes the new value a' k chosen according to the 
equilibrium probability distribution p eq ({cx}). In the case of the Ising chain, the 
transition rates are equivalent to Glauber's ones ^1] . We will use a slightly different 
dynamics consisting in a sequential update of spins. Let us choose a sequence of lattice 
sites G {l,...iV},Vt = nAt,n e IN} and let us define the transition rates in 

discrete time as 

W({a} - {a'},t) = W K{t) ({a} -> {a'}). (18) 

In comparison to Glauber dynamics, only the spin-flip involving the spin o~ k m is possible 
at time t. In the continuous-time limit, the two dynamics are equivalent up to a rescaling 
of time t — > t/N (found for example in the definition of a Monte Carlo step). Indeed, 
when iterating N times the master equation (j3J, one obtains 

p({a}, t + NAt) = (1 - NAt)p({a}, t) 

N-l 

+ At^ P(W}, t) £ W Kit+nAt) ({*'} - {a}) + 0(At 2 ). (19) 

{a'} n=0 

and the Glauber dynamics is recovered if {/c(t + nAt)} n= o...N-i is any circular 
permutation of the set of lattice sites {1, . . . iV}. The equivalence of the two dynamics 
may not hold in the thermodynamic limit iV — > +oo. The time-dependence of 
the transition rates ([18)1 breaks the time-translation invariance of the conditional 
probabilities. However, the effective transition rate in equation (fTH|) is time-independent 
and thus the time-translation invariance is restored in the continuous-time limit if {/«(£)} 
is periodic of period NAt. Again, this may be no more true in the thermodynamic 
limit. In the following, we will assume that {n(t)} satisfies the two above-presented 
conditions, i.e. being periodic of period NAt and that any iV consecutive values are a 
circular permutation of {1, . . . N}. 
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3. Fluctuation-dissipation relation 

3.1. Far-from- equilibrium fluctuation- dissipation relation 

A magnetic field hi is coupled to the spin o~i between the times s and s + At. During 
this interval of time, the transition rates are changed to 



n 

l+k 



(20) 



E CT =±i e 

in order to take into account the additional Zeeman term fihiOi in the Hamiltonian of 
the equilibrium probability distribution (jl3|) . The transition rates are all identical to 
the case hi = apart from the single-spin flip involving the spin <7j. 

Using the Bayes relation and the discrete-time master equation P|). the average of 
the spin at time t > s can be expanded under the following form: 

(a J (t))= 5>;P(M,*) 
M 

= E ^p(W,t|K},s + At)p({or'},s + At) 
= E ^p(W,t|{a'} !S + At)[(l-At)KK},s) 

+ Ai E wJ(.)(K> - K))p(K}^)]- (2i) 
{ CT »} 

W/ 1 being the only quantity depending on the magnetic field in equation (|2*Tj) . only 
remains the second term when «(s) = z after derivating with respect to the magnetic 
field . The derivative leads then to 

~d(a J (t))' 



dh. 



{ CT "> 



dhi 



P(W'},s) 



(22) 



J hi-*0 



This quantity is the magnetisation on site j at time t when an infinitesimal magnetic 
field is applied to the site i between s and s + At, i.e. an integrated response function 
that we will denote Xjiifi I s i s + At]). The derivative of the transition rate Wj 1 defined 
by equation (j2*U|l is easily taken and reads 



dhi 



{{a"} - {a'}) 



hi->0 



0~: 



tanh 



E JikO-' k 

k^i 



(23) 



It turns out to involve the transition rate of the zero- field dynamics ([17)1 . Due to this 
property, the integrated response function can be expressed in terms of thermodynamic 
averages of the zero-field dynamics. Inserting ()23|) into ((22)1 . the integrated response 
function is rewritten as 

Xj i(t,[s;s + NAt})=pAt5 K{sli £ a j p{{a}, t\{a'}, s + At) 
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x 'ai - tanh ((3j2Ji^'k)]Wi({a"} - {a'})p({a"}, s) (24) 

The summation over {a"} can be performed by using the discrete-time master equation 
(J1J). One obtains 

Xji {t,[s;s + NAt}) = (35 K(s)il ]T ^ p({a},t\{a'}, s + At) [a[ - tanh (j3^ W fc ) 

{<x},{<x'} " fc^i 

x p({a'}, s + At) - (1 - At)p({a'}, s)] (25) 

Using a Taylor-expansion of p({a'}, s) in the vicinity of s + At, equation (|25jl can be 
rewritten to lowest order in At as 

Xji (t,[s;s + NAt}) = pAt5 K{s)ti £ ^ t|{<7'}, s + At) [a 1 . - tanh J ife ^ 

p({a'},s + At) + ^({a'},s + At)]. (26) 

The time-translation invariance of conditional probabilities being restored in the 
continuous-time limit, they are function of t — s only and thus satisfy the property 



dP -(W},t\{a'},s) = -^({a},t\{a'},s). 



(27) 



ds dt 
The term involving the time-derivative in equation ()26)1 can thus be rewritten in the 
continuous-time limit as 



P(W},t\W'},s)^({a'},s) 



d_ 

ds 



p({a},t\W},s)p({a'},s) 
-^({a},t\{a'},s)p({a'},s). (28) 



--+^({a},t\{a'}, S ) 



Moreover, the integrated response function Xjii^i [s; s+At]) goes to the response function 
Rji(t, s) in the continuous-time limit : 

rs+At 

Xji(t, [s; s + At]) = J Rji(t, u)du = %(t, s)At + 0(At 2 ). (29) 

Combining equations (J2lj|) . (J28J) and (|2^|). the response function reads in the continuous- 
time limit 



/ d d \ 
R ji (t, s ) = i3(l + - + -Ua j (t) 



(JAS 



Weiss/ x 



k(s),i 



(30) 



where aj^ eiss (s) = tanh ([3J2k^i JikO~k(s)^) is the equilibrium value of the spin <7j in the 
Weiss field created by all other spins at time s. Relation (j3Ti|) generalises equation (HJ). 
The response function Rji(t,s) turns out to be related to time-derivatives of the 
correlation function of the spin o~j at time t with the fluctuations of the spin o~i at time 
s around the equilibrium average af^ eiss (s) of this spin in its Weiss field. In this sense, 
this relation is still a fluctuation-dissipation relation but valid far-from-equilibrium. No 
assumption has been made on the dimension of the space or on the set of exchange 
couplings Jki during the calculation. Moreover, it applies for any initial conditions 
p({a}, 0). The appearance of the prefactor 1 + + is not related to the equilibrium 
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probability distribution of the model but comes from the Markovian properties of the 
dynamics. Generalised response functions are easily calculated along the same lines 
than equation (f3*U|) . The second-order term for example reads 

<(M,r)= ( S^fL ) (t>s>r) (31) 



kji > , h _ u 

( d d\ ( d d d\ 
= P y 1 + ¥s + Ft) r + + ~cTs + di) 6a '^ 6(Ti ^ 6k ^ 6k ^ 

where 5o~j(s) = o~j(s) — crJ Veiss (s). Calculation of non-linear terms requires higher-order 
derivatives of the transition rate as for example 

*8fc-,-)= (*>-) (32) 



Shf s) 



These relations are moreover easily extended to other models. The relations (f3T)|) to 
(132)1 hold for the 0(ri) or the g-state Potts for example where o~i(s) has to be replaced 
by the local order parameter at time s on the site i and <7^ v ' eiss (s) by its average value 
in the Weiss field. Since equations (|3~U|) to ([3*2)1 involve a constraint on the sequence of 
spin-flips, their generalisation to the Ising-Glauber model is not trivial. However, they 
will be of great interest for Monte Carlo simulations. 

3.2. Equilibrium limit 

We will show in this section that the usual expression of the FDT ([IJ is recovered in 
the equilibrium limit. At equilibrium, the probability distribution jp e q({o"}) does not 
depend on time. As a consequence, the integrated response function can be written 
according to equation ("2*""|) as 

X f i {t,[s-s + At}) =/3At5 K{s)tl E 0ip({<7},t|{<7V + At) 

WhW'} 

x 'a[ - tanh J lk a' k )] p cq (K}) (33) 

The hyperbolic tangent can be expressed in terms of the transition ratio of the zero-field 
dynamics (fT7)): 



tanh (/3 g J ik a' k ) p cq ({a'}) = E < -^^-^ 

= E ^({a"} - K}) Pcq (K}). (34) 

Inserting in equation (jHSJ), the integrated response function reads 

{t,[ S] s + At])=(3At5 K{s)>l [ E ^^(M^IM^ + At)^^'}) (35) 

- E ^KW^IKI.s + AtXH/^Kj^K})^^^"}) 



eq 
Xji 
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The first term can be expressed as a thermodynamic average while in the second, one 
needs to get rid first of the transition rate. The Kronecker delta constrains the only 
possible spin-flip to involve site i at time s. As a consequence, Wi({cr"} — > {cr f }) can be 
replaced by W K ( s )({o"} — > {o~'}) and the master equation @ can be applied to equation 
(jHEJ). Moreover, one can show that 

p({a}, t\{a"}, s) = (1 - At)p({a}, t\{a"}, s + At) (36) 
+ At £ W({a"} -> {a'}, s)p({a}, t\{a'}, s + At) + 0(At 2 ). 

This relation is obtained by first putting alone p({a},t\{a"}, s) in the left member 
of the master equation Q and then by iterating the relation to make disappear 
p({a},t\{a"}, s) in the right member. Equation (J36|) is then used to eliminate the 
transition rate from equation (J33j) : 

X$(t,[8-,s + Ab])=l36 lt[a))i [At £ cx, p({a},t\{cx'},s + At)a' lPcq ({a'}) 

W},W'} 

~ Yl ^(M^IKI^ViPcqtK}) (37) 

WhW'} 

+ (1-At) o- 3 p({o-},t\{a'},s + Aty^da'}) . 

The two terms of order At cancel and it remains only 

X %{t, [s; s + At]) = pS^iiatf) [o-tis + At) - ^( s )]) eq (38) 
and in the continuous-time limit, one obtains equilibrium fluctuation-dissipation: 

B%{t,s) = /3^ (s) ,^( ( 7 i (t) ( 7 i ( S )) eq = /3^ W ,i(o- i (t)[tr < (s) - ^ Woiss (s)]) eq (39) 

where the last member is simply equation (j3*Uj) at equilibrium. One recovers the usual 
equilibrium fluctuation-dissipation relation up to a Kronecker delta due the fact that 
the response function is non-zero only for times at which a spin-flip involving the spin 
connected to the magnetic field occurs. 



4. Monte Carlo simulations of the 2d-Ising model 

The discrete-time analogous of expression (jSOj) of the response function enables to study 
the aging displayed by the Ising-Glauber model more accurately than in previous works 
that were based on the numerical estimate of the integrated response function. In the 
first part of this section, the algorithm is given. In the second part, simulations of 
the Glauber dynamics of the two-dimensional Ising model during a quench from the 
paramagnetic phase to the ferromagnetic one are presented. The system is expected 
to display aging, associated with the existence of growing domains corresponding to 
competing ferromagnetic states ^3]. Reversible processes occur in the bulk of domains 
while domain wall rearrangements are irreversible. We will distinguish between quenches 
at the critical temperature T c and below. In both cases, lattice sizes 128 x 128, 256 x 256 
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and 362 x 362 were simulated and the data averaged over 3000, 10000 and 5000 initial 
configurations respectively. For all data, error bars were estimated as the standard 
deviation around the average value. 



4-1. Discrete response function 

During a Monte Carlo simulation, the time is a discrete variable and the time step is 
set to At = 1. Monte Carlo simulations implement indeed the Markov process defined 
by the master equation (j3J) with the choice At — 1. Since simulations are always 
made on finite systems, dynamics with sequential and parallel updates are equivalent 
in the large-time limit up to a time-renormalisation corresponding to the definition of 
a Monte Carlo Step (MCS). The response function can only be defined for continuous 
time processes. However, the integrated response function during one spin-flip o~ t — > a[ 
is the best estimator for the response function Rij(t,s) that we can define. Inserting 
At = 1 into equation (|25j) . the estimator of the response function is simply 

Xji (t, [a; s + 1]) = PS^M® [a^s + 1) - a™ as (s + 1)]) (40) 

where af /eiss (s) = tanh (f3J2k& Jik°~'k{ s ))- m the following, we will be interested only 
on response functions of the form Ru(t,s). In order to reduce statistical fluctuations, 
we have then estimated the response function R(t, s) as the average over all spin-flips 
during one MCS: 

Y N-l 

R(t, s) = — X K (s+n)K( s +n)(t, [s + n; s + n + 1]) (41) 

iV n=0 

The calculation of this quantity is quite simple. Let evolve the simulation until time 
s. For each of the iV next spin-flips cr, — > a' { , store the quantity a[ — af eiss . Note that 
o~\ may be equal to Cj meaning that the system has not changed during this time step. 
However, in strict application of equation (|4()jl . one has nevertheless to store o~i — aj^ eiss . 
After N spin-flips, let the system evolve again until time t. Calculate the response 
function for each site % by multiplying the quantity o\ — af clss stored by the new value 
of the spin <7j and add all these one-site response functions. Repeat the simulation as 
many times as necessary and average the results. The integrated response function can 
be easily calculated by numerical integration of the response function. 

The time-derivative of the correlation function -^Cji(t, s) at time s can be estimated 
by (o~j(t) cjj(s + 1) — (7j(s) ). Again, this quantity is averaged over all spin-flips during 
one MCS. The FD ratio can be estimated as 



R(t, s) 
P£C(t, s) 



Y,n=a{<7 K {s+n)(t) 


^( S+ n)(* +1) "<?„)(* +1) 


"> 


f3Y.n=Q ( a K{s+n){t) 


0~ K (s+n){s + 1) - t7 re ( s+n )(s) 


) 
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Figure 1. Response function of the 2D-Ising model during a quench at the critical 
temperature (inset) and collapse of the scaling function for different values of s. The 
data were obtained with a lattice of size 362 x 362 and averaged over 5000 initial 
configurations. Each curve is surrounded by a clouds of dots corresponding to the 
lower and upper bounds of the error bars. The value a c = 0.115 was used |17| . 



4-2. Quench at the critical temperature 

During a quench at the critical temperature T c , the asymptotic decay of the correlation 
function has been conjectured to be [TBI HI] 

C(t, s) ~ 8—C e (t/ 8 ) (43) 

t,s>l 

where a c = and C c (x) is a scaling function that asymptotically behaves as C c (x) ~^ 
^ j g cr jtj ca j autocorrelation exponent ^H] and z c the dynamical exponent. 
Similarly, the asymptotic behaviour of the response function is 



R(t, s) ^ s- l -^n c {t/s) (44) 

where the scaling function lZ c (x) behaves asymptotically as 1Z c (x) ~ x~ Xc / Zc too. By 

z>i 



integrating over s, one obtains a relation similar to (|44|) for the integrated response 
function that has been checked by large-scale Monte Carlo simulations ^H]- However, 
the relation (}4"4"j) is asymptotic so is not expected to hold for the response function R(t, s) 
with small values of s that are the main contribution to the integrated response function. 
As a consequence, it is difficult to test the asymptotic behaviour of the response function 
in this way. Our approach permits us to avoid these problems and to test directly the 
relation 



The numerical data are plotted in figure ^ For the largest lattice size (L = 362) 
and the smallest value of s (s = 10), errors bars are at most 6 % of the value of the 
response function while for the largest (s = 320), they increase up to 12 %. Indeed, in 
the last case, the response function is of order of l/iV config and so can not be sampled 
accurately. Nevertheless, a fairly good collapse of the data is observed indicating that 
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Cits) 

Figure 2. FD ratio X(t, s) versus t for different values s for the Ising-Glauber model 
quenched at the critical temperature T c . The data were obtained for a lattice 362 x 362 
and averaged over 5000 initial configurations. Each curve is surrounded by a clouds of 
dots corresponding to the lower and upper bounds of the error bars. 



TZ c (t/s) = s 1+a-c R(t, s) is indeed a scaling function (actually we have used (t — s)/s 
instead oit/s but this has no consequence on the asymptotic behaviour). 

We computed the FD ratio X(t, s) using the estimator previously derived and whose 
expression is given by equation (J42)) . The error bars are quite large. The numerical data 
are plotted in figure El for the largest lattice size (L = 362). In contradistinction to 
Cugliandolo conjecture (J2J), the inset of figure El shows that the FD ratio does not 
depend on time only through the correlation function. However, it seems that it may 
be the case in the limit C(t,s) — > 0. On the other hand, it seems that the FD ratio 
depends on time only through t/s and reach a plateau for large enough values of t 
that we may estimate roughly to be ~ 0.33(2). The same value is obtained for 
L = 256 and L = 362 excluding any possibility of finite-size effects. The limit 
has been conjectured to be universal JTj but incompatible values have been given by 
different groups: = 0.26(1) (2H| and X^ = 0.340(5) |23 by Monte Carlo simulations 
and Xoo — 0.35 [221 fo r the 0(l)-model in dimension d = 4 — e. Our estimate is 
compatible with the last two ones. The estimate X^ = 0.26(1) has probably been 
measured for a too-short time t, far from the region where Cugliandolo conjecture ((21) 
and thus equation hold. This puts stress upon the danger of using equation Q to 
compute the FD ratio. 



4-3. Quench below the critical temperature 



The same analysis can be done below T c . In this regime, The correlation function 
decays as FTTj 

C(t, s) tZ>i Mlf{t/s) (45) 
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(t-s)/s 

Figure 3. Response function of the 2D-Ising model during a quench at the 
temperatures T = J/0.6 ~ §T C (inset) and collapse of the scaling function for different 
values of s. The data were obtained for a lattice 256 x 256 and averaged over 10000 
initial configurations. Each curve is surrounded by a clouds of dots corresponding to 
the lower and upper bounds of the error bars. 



where M eq is the equilibrium magnetisation and C(x) a scaling function that 
asymptotically behaves as C(x) ~ x~ x ^ z . The autocorrelation exponent A and the 
dynamical exponent z are expected to take values which are different from those at T c . 
The response function is expected to scale as [TEJ HZj 

R(t, a) t Z>l s-^Uit/s) (46) 

where 1Z(x) ~ x~ x l z . A controversy exists concerning the value of a that has been 
estimated to be either 1/4 or 1/2 j2UI2H]- Our numerical data are presented in the 
figure El We studied lattice sizes only up to L = 256 but calculations were made for 
two temperatures: J/0.6 ~ |T C and J/0.9 ~ T c /2. The error bars are much smaller 
than in the critical case for small values of s. The relative error is at most 2.6 % for 
s = 10 at §T C but increase faster with s: the relative error increases up to 11 % for 
s = 80. As a consequence, the study was limited to the values of s ranging from s = 10 
to s = 80. The response function displays the expected scaling behaviour with 
a = 1/2. However, the collapse is not perfect, especially for the smallest values of s 
but a very small variation of a does not improve it significantly. The value a = 1/4 in 
particular improves the collapse for the small values of s only. The response function 
has probably strong corrections to scaling. Note that corrections have already been 
taken into account for the study of the scaling behaviour of the integrated response 
function [HHH]. 

Combining the relations and (25), the FD ratio X(t,s) = R(t, s)/ (3^C(t, s) 
is predicted to vanish as s~ a below T c . Our numerical estimates for T = J/0.6 ~ |T C 
are plotted in figure |U The statistical errors decrease with the temperature so that the 



Fluctuation-dissip"*""*" — r„„-.,„ m„„.h„„ ui„„ ?„/ 14 



(t-s)/s 

0.001 0.01 0.1 1 10 100 




C(t.s) 

Figure 4. FD ratio X(t, s) versus t for different values s for the Ising-Glauber model 
quenched at the temperature J/0.6/cb — jT c . The data were obtained for a lattice 
256 x 256 and averaged over 10000 initial configurations. Each curve is surrounded by 
a clouds of dots corresponding to the lower and upper bounds of the error bars. 

data are less fluctuating than at T c . As expected, the FD ratio is equal to 1 for small 
values of t — s, signalling that the main contribution to the response function is due to 
equilibrium processes. On the other hand, it vanishes in the limit t ~ s > 1 as s -1 / 2 . 
As shown in figure El the data for s 1//2 X(t, s) collapse for large values of t/s. Moreover, 
figure |U shows unambiguously that the FD ratio does not depend on time only through 
the correlation function. This makes the relation Q invalid. The study of the violation 
of the equilibrium FDT by the usual method relying on the equation (JHJ) would have led 
to erroneous values of the FD ratio. 

5. Conclusion 

Using a formalism similar to Kubo's one in the quantum case, we derive an exact 
expression of the response function of an Ising-Glauber-like model far-from-equilibrium 
(equation (JHUj)). At least for finite systems, the dynamics of our model is equivalent to 
the Glauber dynamics up to a time-renormalisation t — > t/N. The derivation is possible 
because the dynamics consists in a sequential update of the spins and the transition 
rate under a magnetic field can be written as a product of the transition rate without 
magnetic field and of a term depending only on the final spin configuration. The response 
function turns out to be related to time-derivatives of a correlation function involving 
the fluctuations of the spin excited by the magnetic field around its equilibrium average 
in its Weiss field. In this sense, the expression is a generalisation of the equilibrium 
fluctuation-dissipation. Our expression is quite general: no assumption has been made 
during its derivation on the dimension of the space, the set of exchange couplings or 
the initial conditions. Moreover, it can be easily extended to other classical models. 
Generalised and non-linear response functions can be obtained analogously. However, 
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the continuous-time expression (J3U)) may not hold in the thermodynamic limit. Analytic 
results would be desirable. Unfortunately, the response function calculated for the Ising- 
chain by Glauber itself in his original paper does not help because the magnetic field 
was coupled differently to the system (by a multiplicative factor to make the calculation 
feasible while we coupled the field by a modification of the transition rate corresponding 
to the addition of the Zeeman interaction in the equilibrium probability distribution). 
Generalisation to the Ising-Glauber model is not trivial because equation ()39j) sets a 
constrain on the sequence of spin-flips. Nevertheless, It is tempting to imagine that like 
the equilibrium FDT (J3HJ), equations (j3"Uj) to (}3"2"j) hold for the Ising-Glauber model when 
suppressing this constrain on the sequence of spin-flips. 

The expression (|3"Uj) of the response function is then implemented in Monte Carlo 
simulations. Our approach gives access to the response function and the FD ratio 
directly. In particular, the FD ratio can be obtained without assuming the validity 
of the Cugliandolo conjecture (J2J). We then study numerically the homogeneous 
two-dimensional Ising-Glauber model quenched from the paramagnetic phase to the 
ferromagnetic one. Both the response function and the FD ratio display the expected 
scaling behaviour both at T c and below T c . The values, still controversial, of a and 

are estimated to be equal to 1/2 and 0.33(2) respectively, in agreement with 
some previous works. The Cugliandolo conjecture (J2J) does not hold for this model 
apart perhaps at T c in the limit of vanishing correlation functions. This would explain 
discrepancies of previous estimates of relying on Cugliandolo conjecture. The above- 
presented numerical procedure may be extended to many different systems and would 
provided a unambiguous test of Cugliandolo conjecture. We are currently studying the 
dynamics of spin-glasses in this framework. 
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